
global u v w du dv dw a b c dt h m x y dx dy s1 s2 s3 x_best y_best g_best
global Uu Uv Vu Vv Wu Wv

img = imread ('edge-in.png', 'png');
[mu, mv, mw] = size(img);

%m = min(mu,mv);

m = 30;
skip = 2;

u = double(img(1:skip:m*skip,1:skip:m*skip,1));

a = 0
b = -1
c = 0.9
dt = 0.5
h = 1

%m = 30
n = 300
nframes = 300
frameskip = 30
history = zeros(nframes,2);

u0 = zeros(m);
v0 = zeros(m);
w0 = zeros(m);
u = zeros(m);
v = zeros(m);
%w = zeros(m);
du = zeros(m);
dv = zeros(m);
dw = zeros(m);
Uu = zeros(m);
Uv = zeros(m);
Vu = zeros(m);
Vv = zeros(m);
Wu = zeros(m);
Wv = zeros(m);
x = zeros(1,n);
y = zeros(1,n);
z = zeros(1,n);
dx = zeros(1,n);
dy = zeros(1,n);
dz = zeros(1,n);
s1 = zeros(1,n);
s2 = zeros(1,n);
s3 = zeros(1,n);
x_best = 0;
y_best = 0;
g_best = 0;

for i = 1: m
    for j = 1: m
        u0(i,j) = i*h;
        v0(i,j) = j*h;
%        du(i,j) = rand()-0.5;
%        dv(i,j) = rand()-0.5;
    end
end

for k = 1:n
    x(k) = rand()*m*h;
    y(k) = rand()*m*h;
end

x(1) = 30;
y(1) = 20;

%fu(m/2,m/2) = 0.5;
           

close all;
H = figure();
frames = moviein(nframes);

xi = 0;
xj = 0;

for iter = 1:nframes
    
    %u(m/2,m/2) = 5;
    %v(m/2,m/2) = 5;
    %w(10,10) = 5;
    
    for i = 2:m-1
        for j = 2:m-1
            relax(i,j);
        end
    end
    for i = 2:m-1
        for j = 2:m-1
            u(i,j) = u(i,j) + du(i,j)*dt;
            v(i,j) = v(i,j) + dv(i,j)*dt;
            w(i,j) = w(i,j) + dw(i,j)*dt;
        end
    end
    
    [Uv,Uu] = gradient(u, h);
    [Wv,Wu] = gradient(w, h);
    [Vv,Vu] = gradient(v, h);


    
    for k = 1:n
        updateagentedge(k,iter);
    end
    for k = 1:n
        x(k) = x(k) + dx(k)*dt;
        y(k) = y(k) + dy(k)*dt;
        j = round(y(k)/h);
        i = round(x(k)/h);
        if i > 1 && i < m && j > 1 && j < m
            z(k) = u(i,j) + 0.1;
        end
    end
    history(iter,1) = x(1);
    history(iter,2) = y(1);
clf;    
hold on
    %axis([1 m*h 1 m*h]);
    quiver3(x(:),y(:),z(:), dx(:),dy(:),dz(:));
    scatter3(x(:),y(:),z(:));
    %quiver3(u0(:),v0(:),w0(:),Gu(:),Gv(:),Gw(:));
    %surf(w);
    surf(u);
    %plot(history(1:iter,1), history(1:iter, 2));
    axis([1 m*h 1 m*h -5 5]);
    title('edge detector');
    xlabel('u');
    ylabel('v');
    zlabel('w');
    view(-25,80);
hold off

%    subplot(2,2,1);
%    scatter(u(:)+u0(:),v(:)+v0(:),2);%, ceil(1000*ep(:)+100*ek(:) + 1));
%    title('position');
%    axis([1 m*h 1 m*h]);
%    subplot(2,2,2);
%    title('displacement');
%    axis([1 m*h 1 m*h]);
%    subplot(2,2,3);
%   quiver(Gx,Gy);
%    title('gradient');
%    axis([1 m*h 1 m*h]);
    %subplot(2,2,4);
    %contour(Gx.^2.+Gy.^2);
    %title('isoclines');
    %axis([1 m*h 1 m*h]);
   %axis([1 m*h 1 m*h]);
%    frames(:,iter) = getframe(H);
 %   pause(0.1)
    if mod(iter, frameskip) == 0
        frames(:,iter/frameskip) = getframe(H);
    end
end


%clf
movie(H,frames,1,24);
movie2avi(frames, 'out.avi');